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ABSTRACT 

N-body simulations of the hierarchical formation of cosmic structures suffer from the 
problem that the first objects to form always contain just a few particles. Although 
relaxation is not an issue for virialised objects containing millions of particles, coUi- 
sional processes will always dominate within the first structures that collapse. First 
we quantify how the relaxation varies with resolution, softening, and radius within 
isolated equilibrium and non-equilibrium cuspy haloes. We then attempt to determine 
how this numerical effect propagates through a merging hierarchy by measuring the 
local relaxation rates of each particle throughout the hierarchical formation of a dark 
matter halo. The central few percent of the final structures - a region which one might 
naively think is well resolved at the final time since the haloes contains « 10^ particles 
- suffer from high degrees of relaxation. It is not clear how to interpret the effects of 
the accumulated relaxation rate, but we argue that it describes a region within which 
one should be careful about trusting the numerical results. Substructure haloes are 
most affected by relaxation since they contain few particles at a constant energy for 
the entire simulation. We show that relaxation will flatten a cusp in just a few mean 
relaxation times of a halo. We explore the effect of resolution on the degree of relax- 
ation and we find that increasing N slowly reduces the degree of relaxation oc N~^-'^^ 
rather than proportional to N as expected from the coUisionless Boltzmann equation. 
Simulated with the same relative mass resolution (i.e. equal numbers of particles) clus- 
ter mass objects suffer significantly more relaxation than galaxy mass objects since 
they form relatively late and therefore more of the particles spend more time in small 
haloes. 
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1 INTRODUCTION 

A standard technique to study the formation and evolution 
of gravitating systems is to perform an A-body simulation 
in which the mass distribution is discretised into a series 
of softened point particles. This solution can be exact for a 
star cluster where each particle represents a single star, but 
for cosmological simulations of the dark matter each parti- 
cle can be 10^" times larger than the GeV mass candidates 
being simulated. In this approach the particles represent a 
coarse grained sampling of phase space which sets a mass 
and spatial resolution. Unfortunately these super-massive 
particles will undergo two body encounters that lead to en- 
ergy transfer as the system tends towards equipartition. In 
the real Universe the dark matter particles are essentially 
coUisionless and pass unperturbed past each other. 

The processes of relaxation is difficult to quantify, but 
in the large N limit the discreteness effects inherent to the 
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N-body technique vanish, so one tries to use as large a num- 
ber of particles as computationally possible. Increasing the 
mass resolution of a given simulation allows a convergence 
tes t of propertie s such as the dark matter de nsity profile 
i.e. iMoore et alj iITqO S'). Ghigna ct al. (2000), K lvpin et all 
J200lf) and lPower et alJ ||2003D . These authors find that to 
resolve the central one per cent of a dark matter halo the 
entire system must contain of the order a million particles. 
It is not known what process sets this resolution scale since 
with one million particles relaxation is expected to be small, 
even at one per cent of the virial radius. 

Unfortunately in most cosmological simulations the im- 
portance of two body interactions does not vanish if one 
increases N . Structure formation in the cold dark matter 
(CDM) model occurs hierarchically since there is power 
on all scales, so the first objects that form in a simula - 
tion always contain only a few particles jMoore et alj|200lll . 
(jBinnov fc Kneba .2002') . 'With higher resolution the first 
structures form earlier and have higher physical densities 
because they condense out of a denser environment. Two 
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body relaxation increases with density, so it is not clear if 
increasing the resolution can diminish the overall amount of 
two body relaxation in a CDM simulation, i.e. if testing for 
convergence by increasing the mass resolution is appropri- 
ate. 

In isolated equilibrium systems relaxation rates can be 
meeisured from the energy dispersion. In cosmological sim- 
ulations one can measure the amou nt of mass segr egation 
of multi-mass simulations ( Binncv fc Kneb3 | 2.002j) where 
lighter particles gain more energy from coUisional processes 
than the heavier particles. In Section|5|we present a Fokker- 
Planck type relaxation time estimate, which was fitted to 
a series of test simulations (Section 13.111 where we explore 
the relaxation rates as a function of N, radius and softening 
parameter in both equilibrium and non-equilibrium cuspy 
haloes. We then use this local relaxation rate estimate to 
follow the relaxation history of each particle during 1000 
time-steps of a hierarchical CDM simulation. The resulting 
degree of relaxation as a function of spatial position within 
galaxy and cluster mass haloes is analysed in Section 0] and 
in Section |3 we discuss the effects that relaxation has on 
haloes at 2 = 0. 



2 A LOCAL RELAXATION TIME ESTIMATE 

In thi s paper we adopt the energy definition of the relaxation 
time l|Chandrasekharll942h stating that the mean relaxation 
time T of a group of stars is the time after which the mean 
square energy change due to successive encounters equals 
the mean kinetic energy of the group: 
-2 
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where AE{At) is the energy difference of one particle after 
time interval At and the bar denotes the group average. 

Note that this time is half of t he relaxation time r„ 
defined in lBinnev fc Tremaind (j^S^) who calculate a mean 
velocity change, because AE^/E^ ~ 2Av^/v^. The orbit 
averaged Fokker-Planck estimate for r„ is 



T„ = 0.34- 



(2) 
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where a is the one dimensional velocity dispersion, p the 
density and m the particle mass. The parameters bmin and 
bmax are the minimum and maximum limits for the impact 
parameter. 

To assess the degree of relaxation in cosmological sim- 
ulations (section |1J we estimate the local relaxation rate for 
each particle after every time-step and integrate this up over 
the whole run: 

k 

d{tk) ■-^rLE{t„)At. (3) 

n = l 

For the local relaxation rate we use a formula similar to ^ 

TLE = — ^l ^f ,7 = 0.17. (4) 
rLE G^mpC 

The value of 7, and the parameters in the Coulomb loga- 
rithm C are chosen to roughly fit measured relaxation times 
of equilibrium haloes, see section The Coulomb loga- 
rithm is 



ln(l + A^ 



A^ 
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the analytical calculation for Newtonian potentials shows 
that bmin = bo — 2Gm/v'^^i, bo is the impact parameter 
where the defiection angle reaches 7r/2 ijBertinl lioOO'l. In a 
softened potential the scattering calculation has to be done 
numerically and the results agree roughly with the Newto- 
nian case if one sets bmin ~ e, i.e. one ignores all encounters 
with an imp act parameter smaller than the softening length 
jTheijll998ll . We set 



max(Cm/3a^, e) ~ max(bo, e) , 



(6) 



because ii^^j = 6a^. The proper choice of bmax is controver- 
sial, it is not clear whether it should be related to the size 
of the whole system or to the mean interparticle distance. 
For cosmological simulations we prefer the second choice 



/3(m/p)^/^ 



(7) 



because this is a local quantity that is easy to measure and 
less ambiguous than defining the size of irregular shaped, 
collapsing structures. 

We calculate the local velocity dispersion and density 
surrounding each particle by averaging over its 16 nearest 
neighbours. We do a simple top-hat average, because us- 
ing an SPH spline kernel leads to biased results when using 
only 16 particles. We found good agreement with all mea- 
sured relaxation rates when using /3 = 10 and 7 = 0.17. 
Averaging over different numbers of nearest neighbours the 
optimal parameters differ slightly due to different amounts 
of numerical noise in the local density and velocity disper- 



3 TWO BODY RELAXATION IN SPHERICAL 
HALO MODELS 

In this section we present a number of test cases which we 
used to gauge our local estimate Q for the degree of relax- 
ation and show that it agrees quite well (within 15 per cent) 
with measured levels of relaxation for a wide range of parti- 
cle numbers, softening lengths and virial ratios. This range 
covers the haloes that form in cosmological simulations, 
however all the test cases are isotropic, spherical haloes. 
Haloes in cosmological simulations are close to isotropic, but 
are triaxial and contain substructures. But one can argue 
that locally the two are similar, and if a local relaxation 
time estimate works in the spherical haloes, it should give a 
reasonable estimate also in the cosmological case. 



3.1 Equilibrium haloes 

In an equilibrium model the energy of each particle would be 
conserved in the A'' — > cxd limit. For finite A*' the energies of 
particles suffer abrupt changes due to encounters. Therefore 
we just have to measure these energy differences to get the 
relaxation time with Q. 

Here we present a sequ ence of tests usi ng spherical and 
isotropic Hernquist models jHernauistll990h which are a rea- 
sonable approximation to the haloes found in cosmological 
N-body simulations. (We f ound no difference between these 
simulations and tests using iNavarro. Frenk fc White! (Il996l) 
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and lMoore et all ilQQgl) profiles, constructed by solving for 
the exact p hase space distribution function numeric ally as 
described in lKazantzidis. Magorrian fc Moorel (|200.'^) .'l The 
density profile of the Hernquist model is 



We set ^ the total mass to M = 3.5 x lO^M© and the scale 
length to a = 10 kpc. Then the half mass radius is — 
2.4a = 24 kpc and the crossing time at half mass radius is 
Tc = rh/vcirc(rh) — 1.3 Gyr. 

All the simulations have been carried out using PKD- 
GRAV, a state of the art, multi-stepping, parallel tree-code 
ijStadel 2nOLy The time-steps are chosen proportional to the 
square root of the softe ning length over the acceleration on 
each particle, Ati = "q^ftfai. We use x] = 0.25, and a node- 
opening angle Q — 0.55 for all runs in this section, expect 
the long term integrations in subsection 13.31 where we use 
rj = 0.03. Energy conservation was better than 0.1 per cent 
after several crossing times for all runs in this section. 

Due to softening the initial models are not exactly in 
equilibrium, the total kinetic energy is a few percent larger 
than half of the potential energy. For this reason we evolved 
the models for five crossing times before measuring the en- 
ergy dispersion, which results in up to 10 per cent longer 
relaxation times. 

Figure shows AE'^{At) as a function of time within 
haloes constructed with A'^ — 10* and N = 100 particles. 
The upper panel shows that for small numbers of particles 
N < 10^ AE{At) becomes very noisy since there are fewer, 
but more significant encounters. To obtain more reliable re- 
sults in small N groups we added a sufficiently large number 
(10*) of massless tracer particles following the same distri- 
bution in real and velocity space as the massive particles. 
The open squares in Figure Q show the 'energy' dispersion 
of the tracers, which in large N groups is just the same as 
the energy dispersion of the massive particles, but it evolves 
much smoother with time in small N groups. We obtain the 
mean relaxation times Q of these haloes with linear fits to 
the energy dispersion of the tracer particles (open squares), 
taking into account points where AE'^ < 0.2, i.e. in the up- 
per panel only the points left of the vertical bar, to make 
sure that At is small compared to the relaxation time. The 
local relaxation estimate (solid line) (jlj gives similar aver- 
age degrees of relaxation in these test cases (see also Figures 
HandlHJ. 

Note that the tracers are not in equilibrium with the 
halo, on average they gain speed in encounters and are 
ejected from the core. Typically after one relaxation time 
the number of tracers inside of r = a drops to one half of 
the initial number. Therefore it is important to use a At 
shorter than T and to add the tracers after evolving the 
halo for five crossing times, otherwise the relaxation in the 
core is not sufficiently refiected in the energy dispersion of 
the tracers and T could be underestimated significantly. 



^ To rescale the results to different size haloes just change the 
distance scale by some factor x — > fx, mass scale M 
and the dynamical and relaxation timescales do not change. To 
rescale to different timescales T — > cT do M — > c~^M with fixed 
length scale. 
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Figure 1. Mean squared energy change AE^{At) as a function 
of time. The filled squares show the mean squared energy change 
of the actual particles, the open squares for the massless tracer 
particles. The dotted lines are linear fits to the mean squared 
energy change of the tracers, in the upper panel only the points 
left of the vertical bar are taken into account. The solid lines are 
averages of over all massive particles. 

3.1.1 Dependence on N and e 

Figures 121 and |5] show the measured mean relaxation times 
as a function of TV and softening parameter, e, compared 
with the average over all particles of the local relaxation 
estimate I^J- We find that the measured relaxation times 
are proportional to N (dashed reference line), rather than 
to N/\n{N) (dotted line), as expected for a softened po- 
tential. The same result wa s foun d for King models by 
iHuang. Dubinski fc Carlberg (199^. The local estimate of 
the relaxation time increases slightly faster with N than the 
measured values. This is due to the fact that we choose a 
maximum impact parameter proportional to the mean in- 
terparticle separation, i.e. bmax oc N~^^^ , but the difference 
is less than 10 per cent for all relaxation times shorter than 
a Hubble time, i.e. for N ^ 10*. 

The dependence on the softening length is shown in 
Figure 13 for a N = 10* model. The measured values (filled 
squares) increase faster with e than the loca l estim ate (open 
squares ) (like in Figure 2 of iHuang et"al] Il993ll . but the 
differences are small (^15 per cent) for realistic softenings 
e<^0.1a — 1 kpc. The average relaxation time of this model 
increases from 30 Gyr to 180 Gyr when the softening param- 
eter is changed from 0.01 kpc to 1 kpc. This is slightly higher 
than expected from the scaling with ln(e) since the density 
profile and central cusp are better resolved with smaller soft- 
ening. 



3.1.2 Dependence on radius 

Measuring the relaxation time as a function of radius T{r) 
proved to be quite difficult. The most credible method seems 
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Figure 2. Average relaxation times of isotropic Hernquist models 
versus particle number A'^, with a constant softening e = 0.5 kpc. 
The filled squares are the measured relaxation times, v^ith error 
bars form the linear fit of AE'^{At). The open squares are the 
local estimates of the relaxation time El . 
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Figure 3. Average relaxation times of an Af = 10* Hernquist 
models versus softening length e. The filled squares are the mea- 
sured relaxation times, the open squares are the local estimate 



to take Q and replace the average over all particles by the 
average over those which are in the corresponding radial 
bin at the beginning (or at the end) of the time interval At. 
Clearly one has to choose At <^ 5r/a{r), where 5r is the size 
of the bins, to make sure that most particles spend most of 
At in the same bin. This could also be achieved by placing 
the tracers on circular orbits, which leads to very similar 



results for T(r) if r ru, but in to the centre this method 
fails, because there the circular velocity is much smaller than 
the velocity dispersion^. 

In Figure @ we plot the relaxation rate against radius 
for a Hernquist model with lO** particles. We measured the 
energy dispersion (filled squares) during At — 0.1 Gyr for 
each particle, and averaged the values of particles starting in 
the same radial bin. We also measured the local relaxation 
rate tle © at 100 time-steps during At for each particle 
and summed them up. The radial averages are plotted with 
open squares. Again the agreement with the measured en- 
ergy changes is better than 35 per cent except in the last 
three bins where the relaxation rates are many thousands 
of Gyrs and the local relaxation measurement overestimates 
the true rate of relaxation. In the inner three bins the cross- 
ing times are shorter than 0.2 Gyr, i.e. many of the particles 
had time to move through these bins during At = 0.1 Gyr. 

The dashed line is the inverse of half the Fokker-Planck 
estimate ^ with A = r^/e, calculated using all particles 
in the bin, not only from 16 nearest neighbours. It scales 
like the phase space density p{r)/a^{r) which scales almost 
exactly like in a Hernquist model. The local estimate 
also follows this scaling in the seven outer bins. The 
measured relaxation scales more like oc r~'^ in the outer re- 
gion, but the slope depends strongly on At, i.e. on how many 
particles from the core with 100 times higher relaxation rate 
have had time to reach the outer region. 

The average relaxation times^ are about 10 times 
shorter than measured relaxation times at half mass radius, 
due to the fast relaxation in the high density core of cuspy 
haloes. For a less concentrated King model (^o = 5) the 
Fokker-Planck estimate seems to agree not only with the 
measured relaxation time at half mass radius , but also with 
the mean relaxation time jHuang et alJll993^ ■ 



3.2 Non-equilibrium systems 

The definition of the relaxation time Q is mostly used for 
systems close to dynamical equilibrium like globular clus- 
ters, i.e. systems with constant (or only slowly changing) 
mean kinetic energy. In this case the mean relaxation rate is 
also (roughly) constant and the degree of relaxation grows 
linear with time (like in Figure 0. Generalising this defi- 
nition to non-equilibrium situations is straightforward: The 
degree of relaxation of a group at some given time is the 
accumulated mean square energy change due to encounters 
divided by the mean square kinetic energy at this time. 
Now not all energy changes are due to encounters, so 



^ From a convolution of the velocity distributions one finds that 
the relative velocities in encounters with tracers on circular orbits 
have a different dispersion (tr^ ^ = 3cr^ + v'^^^^) than those of 
encounters between the massive particles (ct^ ^ = 6ct^). In non 
isothermal haloes, ^ ^circ ^^'^ ''^'^ leads to systematic errors 
in the measured relaxation times. In Hernquist models 3(J^ ~ 
^circ holds only for r r/^, and indeed we found good agreement 
with the other methods only in this range. 

Note that analytically the average of the relaxation rate esti- 
mate r^E divergent for models with central cusps oc 
and steeper: T = ri^^{r)p{r)T^ dr oc p{r)r^ dr = 
/(f p{r)dr. 



© 2003 RAS, MNRAS OOO.ITHTTI 



Two body relaxation in CDM simulations 5 



-1 1 1 1 1 1 1 M 


\ 














. " \ 








\- 

: °\ 






- 


_ T,, = 105 Gyr 


\ 




- 


- T = 1 13 Gyr 




v 


- 






\ 

■ \ □ 

\ 








\ 

■ 


\ ° i 

\ - 


N = 10* 








_ E = 0.5 kpc 






■ : 



1 10 

[kpc] 

Figure 4. Relaxation rate of an Af = 10^ halo vs. radius. The 
filled squares are the measured relaxation, the open squares arc 
the local estimate r^^g, calculated from 16 nearest neighbours 
during 0.1 Gyr. The horizontal lines give the halo averages of 
measured (solid line) and estimated (dotted line) relaxation rates. 
We also plot 2/Tv (dashed line) calculated from all particles in 
the radial bin. 
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Figure 5. Relaxation times of a TV = 10* non-equilibrium Hern- 
quist model vs. virial ratio. The filled squares are measured with 
/ tJ^ , the open squares are the local estimate J4} . The reference 
line is oc a^'^ oc ifl I p. 



one needs another method to measure the degree of relax- 
ation. Again we use massless tracer particles like in the last 
section. Instead of setting up their initial conditions exactly 
like those of the massive particles, one can also restrict the 
tracers to a common orbital plane. In a spherical system this 
does not change their density profile nor the relative veloc- 
ities in encounters with the massive particles. If we choose 
the orbital plane of the tracers to be the a;j/-plane, then the 
accumulated energy change due to encounters Aenci? is 

Aencg(f)^ _ Acncgfc.n(t)^' _ ^ Aenc-»(t)^ ^ ^ v{t)l 

as long as Uz(t)^ <^ i-e. for small degrees of relaxation. 

This relates the energy dispersion to a more demonstrative 
quantity and in the edge on view of the xy-plane one can 
actually observe the relaxation process since the degree of re- 
laxation is roughly proportional to the thickness of the disk*. 
When the amount of relaxation approaches unity tends 
to underestimate the degree of relaxation, because tracers on 
new out of plane orbits will eventually reach a turnaround 
point where Vz = Q. Also the probability that one tracer 
suffers more than one encounter grows with time, therefore 
one should place a new set of tracers into the plane to get an 
accurate relaxation rate as soon as the amount of relaxation 
is close to unity. 

With ® we can measure the amount of relaxation 
in non-equilibrium situations, we only need one symmetry 
plane to be able to apply this method, therefore it allows us 

* Simulation movies are available at: 

www-theorie.physik.unizh.ch / ~dicmand / tbr / 



to measure the amount of relaxation during a collapse or a 
merger. 



3.2.1 Collapsing haloes 

In CDM simulations the virial ratio a = 2Ekin/\Epot\ is 
close to zero at the beginning of a halo collapse and grows 
towards unity as the halo reaches dynamical equilibrium. In 
the previous sections we showed that the local estimate @ 
works for a = 1, but for q — > the phase space density 
oc a~^'^ goes to infinity. Down to which virial ratio can we 
trust our local estimate? 

To answer this question we began with equilibrium 
Hernquist models (same parameters as in section 15.1. in with 
10* particles and multiplied all velocities with y/a. There- 
fore oc a^'^ and the phase space density p/a^ oc q~^'^ . 
To these models we added 10'' massless tracer particles with 
the same phase-space distribution and tilted their orbital 
planes into the sy-plane. Then we measure vljlP" and the 
local estimate l|3J at 100 time-steps during the first 0.1 Gyr 
of the collapse. Linear fits give the relaxation times plotted 
in Figure the dashed line shows the scaling oc a*'^ ex- 
pected from the Fokker-Planck type estimates and @, 
since these times are proportional to one over phase-space 
density. Our relaxation time estimate becomes very small 
when the virial ratio goes to zero, but the measured relax- 
ation times remain on the order of a few dynamical times. 
Tle is within a factor of two for a>0.075 and within 10 per 
cent for a<;0.5. Also after the first 0.1 Gyr the local esti- 
mate follows the measured values, one example (a = 0.25) 
is plotted in the top panel of Figure 
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Figure 6. Average degrees of relaxation during a collapse of an 
N = 10* Hernquist halo, the initial virial ratio a is 0.25 (top). 
The lines are the average of the local relaxation estimate , the 
squares are the amount of relaxation measured from a group of 
tracer set in a plane initially j5J . The bottom panel shows average 
degrees of relaxation during a a merger of an Af = 10"^ halo into 
a ten times more massive system with N = 10*. Here the open 
squares are the measured average for the satellite and the filled 
squares for the main halo. 



3.2.2 Mergers 

The last test case for the local relaxation estimate is a 
merger of a small system into a more massive one. The prob- 
lem in this case is that the small halo gains a lot of kinetic 
energy when falling into the main halo, so its accumulated 
energy changes due to encounters can become smaller rela- 
tive to the mean kinetic energy of the group. A local estimate 
can never capture this decrease since it can not know about 
the gain in external kinetic energy. For a surviving subhalo 
one can argue that its accumulated energy changes should 
still be compared to its roughly constant internal kinetic en- 
ergy to get an estimate of how affected it is by relaxation. 
But for haloes that are disrupted (and stripped particles 
from subhaloes) one has to worry about how their overesti- 
mated degrees of relaxation affect the average relaxation of 
the main halo. 

The bottom panel of Figure (|SJ shows how relaxation 
develops in a head on merger of a. N = 10^, a = 3 kpc halo 
into a ten times more massive halo with A'^ = 10* and a — 10 
kpc. The initial separation was 100 kpc with a small relative 
velocity of 2.2 km/s. The satellite falls in and reaches the 
centre of the main halo after 3 Gyr. Note how the relaxation 
of the satellite (open squares) decreases during infall, this 
can not be followed by the average local estimate I^J of the 
satellite (dashed line), but still the average over the whole 
system (solid line) gives a good estimate for the mean degree 
of relaxation. We also verified this for equal mass mergers, 
there the decrease during infall is small because both haloes 
have quite large internal energies initially. 
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Figure 7. From top to bottom, the mass fraction within 
0.89, 1.58, 2.81, 5.00, 8.90, 15.8, 28.1 and 50.0 kpc of a Af = 4'000 
Hernquist models vs. time. 



3.3 Evolution of isolated Halos 

The dynamical evolution of globular clusters is driven by 
relaxation, which can lead to core collapse and evaporation 
on a timescale of a few tens of half mass relaxation times, 
i.e. the core loses energy to an e xpanding outer envelope of 
stars and get s denser and hotter (iBinnev fc TremaindllQs'^ : 
ISpitzeillT987i) . 

In the next section we show that haloes in cosmologi- 
cal simulations are typically between one and ten mean re- 
laxation times old. (In terms of the much longer half mass 
relaxation time they are younger than one or two half mass 
relaxation times.) Therefore we do not expect that density 
profiles in cosmological simulations are significantly affected 
by the core collapse process. 

Studies of globular cluster evolution start with mod- 
els that are isothermal in the centre (e.g. Plummer spheres, 
King models) and then show a slow but monotonic density 
increase in the core. In contrast the cuspy haloes in cosmo- 
logical simulations are not isothermal: the velocity disper- 
sion decreases in the central regions. In this case relaxation 
leads to an energy flow inwards, the core evolves towards a 
less dense, isothermal state first. Later the system evolves 
just l ike the models in globular cl uster calculatio ns ( Quinla3 
Il99d) . The N-body simulations of lHavashi et al.l EoOSi] show 
this evolution starting from an NFW profile. We confirmed 
their result by evolving an A'^ = 4'000 Hernquist model for 
360 crossing times (see Figure |7J. 

Figure |H| shows the evolution of five A'^ = 4' 000 Hern- 
quist models during 16 Gyr, using a softening of e = O.lkpc. 
For this long term evolution we use a more conservative time 
step parameter r] = 0.03, energy is conserved within 0.36 per 
cent even after 360 crossing times. The crossing time at the 
half mass radius is 1.3 Gyr, the initial mean relaxation time 
is 16 Gyr and the initial half mass relaxation time is 71 
Gyr. After 50 Gyr this halo is about three mean relaxation 
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Figure 8. From top to bottom, the mass fraction within 
0.89,1.58,2.81,5.00,8.90,15.8,28.1 and 50.0 kpc of five N = 
4'000 Hernquist models vs. time. 
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Figure 9. Averaged density profiles of five TV = 4'000 Hernquist 
models (same models as in Figure ISl. initial profile (dashed line), 
after 16 Gyr (dotted line) and after 50 Gyr (solid line). The points 
indicate the radius of the outer borders of the spherical bins. 



times old, a realistic value for haloes in current cosmological 
simulations, see section (|1J. The same would happen to a 
N = 100 halo in only 1.25 Gyr, we use iV = 4'000 just to 
have a well defined density profile down to 0.1a — 1 kpc. 

Within few mean relaxation times the velocity disper- 
sion rises in the centre and at the scale radius it drops 
slightly, in the end the core is already close to isothermal. 
This energy gain is compensated with core expansion, the 
inner two mass shells plotted in Figure |H| clearly lose mass. 



This is not a numerical artifact, the softening we used is 
much smaller than the inner bin and in a iV = 4 X 10'' refer- 
ence model the mass loss in the inner two bins is about ten 
times slower, i.e. this is really an effect driven by relaxation. 
The corresponding density profiles are less steep in the in- 
ner 3 per cent of the halo, see Figure |U] and show constant 
density cores in this region. 



4 RELAXATION IN COSMOLOGICAL 
SIMULATIONS 

Here we present results from four low to medium resolu- 
tion ACDM simulations (^a = 0.7, fi™ = 0.3, as = 1.0). 
We generate initial conditions with the GRAFIC2 package 
jBertschingeill200ll) . We start with a 128^ particle cubic 
grid with a comoving cube size of 60Mpc (particle mass 
3.6 X IO'^Mq). Later we refined two interesting re- 
gions, in the first one a cluster halo (M200 = 7.4 x IO^^Mq, 
'"200 = 1440 kpc) forms, the refinement factors are 2 and 3 in 
length, ie. 8 and 27 in mass (run C2,C3). In the second region 
a galaxy size halo (M200 ~ 1.4 x lO^^M©, r2oo = 350 kpc) 
forms. There we used a refinement of 9 in length, ie. 729 in 
mass, and included a buffer region, about 2 Mpc deep, with 
an intermediate refinement factor of 3 (run G9). We start 
the simulations when the standard deviation of the density 
fiuctuations in the refined region reaches 0.2. The softenings 
used in the refined regions are e — 1.86 kpc for the cluster 
and e = 0.5 kpc for the galaxy, i.e. e ~ 0.0013ri,ir in both 
cases. We also run the unrefined cube again with this soft- 
ening in the cluster forming region (run CO). The numerical 
parameters are as in section ITTl but at late epochs we use a 
larger node-opening angle to speedup the runs, 9 = 0.7 for 
z <2. 

In non equilibrium, non spherical haloes of a cos- 
mological simulation it is not possible to measure two 
b ody relaxation times wi th the methods used in section 
IHl iBinnev fc Knebd ll2002l) used initial conditions with two 
species of particles with different mass. Both species start 
from a regular lattice, such that the nodes of one grid are 
at the centres of the cells of the other, and both are then 
displaced according to the Zel'dovich approximation. In a 
coUisionless simulation the final distribution of the particles 
would be independent of mass. They found differences in 
the number density of the light and heavy particles in the 
centres of haloes. 

Here we compliment the study of Binney & Knebe by 
applying the results of the previous sections to cosmological 
simulations. We assign a degree of relaxation d to each par- 
ticle, which is calculated after each of 1000 fixed time-steps 
At from the local relaxation rate estimate @ and summed 
up over the whole cosmological simulation Q. 

As shown in section 1!). 2. II tle reflects the measured re- 
laxation rates only for virial ratios a^O.l. Since this is not 
case for the first steps in a CDM simulation, we set tle to 
zero before the local density reaches some threshold. When 
the local density reaches 6po (— density at turnaround in 
the spherical collapse) the typical values for a are close to 
0.4, later at 170po (— density at virialisation in the spheri- 
cal collapse) a is close to unity. We used these two density 
thresholds, in the first case we write dxA for the "degree of 
relaxation since turnaround", otherwise dvm for "degree of 
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Table 1. Average Dej 


;rees of Relaxation 






Run 


CO 


C2 


C3 


G9 




20' 500 


177'000 


650'000 


250'000 


do inside 0.1r200 


5.98 


3.62 


3.06 


1.67 


UQ iiitiiut; / 200 


5.23 


3.34 


2.52 


1.15 


dxA inside 0.1r200 


4.74 


2.40 


1.78 


0.72 


dxA inside r2oo 


3.67 


2.42 


2.12 


1.02 


dvm inside 0.1r200 


3.58 


1.61 


1.17 


0.42 


dviR inside r2oo 


2.50 


1.52 


1.34 


0.58 



relaxation since virialisation" . The relaxation averages over 
all particles inside r2oo and 0.1r-2oo at 2 = are given in 
table El 

4.1 Number of particles 

A reassuring result is that with more particles the simula- 
tions are less affected by two-body relaxation, even though 
one resolves more small A'^ progenitors. This confirms the 
significance of convergence tests that vary the number of 
particles. The average degree of relaxation inside of 0.1r2oo 
scales like N~'^'^ , and the relaxation inside of 7-200 hke A*""'^. 

Figure [TTI shows the relaxation in the cluster as a func- 
tion of the final particles position, for three different resolu- 
tions. In the outer part (r>0.1r2oo) the cluster has substruc- 
ture, which are small A'^ systems that exist at the present 
time, so they are still relaxing at a high rate at 2: = (bot- 
tom panel). Substructure haloes with A'' ~ 500 can reach 
averages of dym — 10 in all runs, the highest peaks in dyiR 
are found in the centre of substructure haloes, where dym 
can be as high as 100, much higher than in the centre of the 
host halo (see Figure [TUl . 

Note that the degrees of relaxation (top and middle 
panel) are much larger than what you would estimate us- 
ing the final distribution of particles (bottom panel) . Other 
studies consider only the relaxation rate at z = 0, and claim 
to resolve a halo down to a radius wher e this relaxation ti me 
^Lisi^ = 0) is larger than Hubbl e time jPower et alJ2003ll o r 
larger than three Hubble times llFukushig ^^MakirMl200ll) . 
This radius scales oc A'^""'^, where as convergenc e in N- body 
simulations seem to be slower; In iMoore et alj lll998h and 
iGhigna et alJ i2000() the resolved radii are determined by 
comparing density profiles between simulations with differ- 
ent numbers of particles. They found that the "resolved 
radius" r ~ 0.5(A'^20o/V2oo)~^''^ for a wide range of A^200 
fromm 10^ to 10^ ^ It appears like the resolved radius scales 
in the same way as the average degree of relaxation, but 
further relaxation studies for a wider range A'' are needed to 
verify this. 

4.2 Mass and time dependence 

Figure fT?l compares the relaxation rate within the high res- 
olution cluster and the galaxy simulations. The average de- 
gree of relaxation at 2 = for the galaxy (dyiR — 0.58) 
is much smaller than for the cluster (dyiR ~ 1.34), even 
though A'^ is larger for the cluster and therefore the present 
relaxation rate tle^z = 0) is smaller in the cluster. The 




Figure 10. Maps of the cluster's density (top) and relaxation 
since virialisation (bottom) out to r200 for the C3 run {N200 — 
650'000). The logarithmic scale for the degree of relaxation goes 
form 0.01 (black) to 100 (white). 



reason is that the cluster forms much later than the galaxy 
therefore most of its particles have spent a longer period of 
time in small N progenitor haloes. 

In Figures [T^ and n~T1 we plot how the degree of relax- 
ation increases with time for both haloes for particles within 
10 per cent of the virial radius and for all the particles within 
the virial radius. Most of the relaxation within the central 
region of the galaxy occurs within the first couple of Gyrs 
of the evolution of the Universe. The cluster forms over a 
longer timescale and this is refiected in the longer increase 
in the degree of relaxation with time. 

The cluster runs (C0,C2,C3) show how relaxation in 
small A' groups starts earlier, after 1 Gyr the highest reso- 
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Figure 11. Relaxation versus radius. The solid line is for run C3, 
the dashed line for run C2 and the dotted line for run CO. 





Figure 12. Relaxation vs. radius. The solid line is the cluster 
run C3 (N^i,- ~ 650'000), the dot - dashed line for the galaxy 
(N^ir ~ 250'000). 



lution run (C3) is most affected by relaxation. Tliis result is 
not an artifact from using a density threshold, we checked 
that the do shows the same behaviour as dviR- The en- 
tire haloes (bottom panel) show some relaxation during the 
whole simulation which arises from the poorly resolved sub- 
structure haloes in the outer regions. 



5 CONCLUSIONS AND DISCUSSION 

A'^-body cosmological simulations attempt to model a col- 
lisionless system of particles using a technique that is in- 
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Figure 13. The degree of relaxation, rfyiRi averaged over all 
particles as a function of time. In the top panel we average over 
particles inside 0.1 r2oo at 2: = and in the bottom panel over all 
inside r200- 
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Figure 14. Same as Figure ITsl but as a function of rcdshift. zq 
is the starting rcdshift of the runs. 



herently coUisional on small scales. We have examined the 
relaxation rates of isolated equilibrium cuspy haloes as a 
function of particle number, radius and softening parame- 
ter. Our results apply primarily to n-body, P^M codes, such 
as direct or treecodes. Adaptive grid based methods, such 
as ART jKravtsov. Klvpin fc KhokhlovUlQQTTl and MLAPM 
jKnebe. Green fc Binnevll200 j) also seem to suffer from re- 
laxation but at slightly lower rat es than from those quanti- 
fied here (iBinnev fc Knebdl2002|) . 

We show how one can define a local relaxation timescale 
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for each particle by measuring its local phase space density 
which we applied to cosmological simulations in an attempt 
to determine the regions most affected by numerical relax- 
ation. We summarise our results here: 

(i) The relaxation rates in cuspy dark matter haloes are 
in good agreement with the rates predicted by the orbit av- 
eraged Fokker-Planck equation. However the average relax- 
ation time is an order of magnitude less than that measured 
at the half mass radius. 

(ii) We verify that the average relaxation time of a halo 
is proportional to the number of particles it contains and to 
the inverse natural logarithm of the softening parameter. 

(iii) The relaxation time is proportional to the local phase 
space density which allows us to measure the cumulative 
amount of relaxation each particle undergoes during the evo- 
lution of a halo. 

(iv) We show that we can measure the relaxation rate 
in collapsing or non-equilibrium haloes that have kinetic to 
potential energy ratios up to ten times smaller than the 
equilibrium value. 

(v) Averaging over several simulations of cuspy Ifernquist 
haloes we show that within few mean relaxation times the 
central cusp is transformed into a constant density core. 

(vi) We show that the hierarchical build up of galaxy 
or cluster mass haloes leads to a greatly enhanced degree 
of relaxation within their central regions. The substruc- 
ture haloes suffer from the highest rates of relaxation since 
they contain the fewest particles for the longest period of 
time. Subhaloes are typically several relaxation times old 
therefore one should be cautious about interpreting their 
internal structure using simulations of order lO'^ particles 
JStoehr et. aJl2002h . 

(vii) Cluster haloes suffer three times the amount of re- 
laxation as galaxy haloes simulated at the same relative spa- 
tial and mass resolution. This is because the cluster forms 
later and more of its particles spend time in poorly resolved 
progenitors. 

(viii) Increasing the resolution (A'') at a fixed force soften- 
ing reduces the the accumulated amount of relaxation. The 
average degree of relaxation in CDM haloes at z = scales 
oc N-°-'^, in the inner 10 percent oc A'' ^'^ . Relaxation may 
therefore provide a simple explanation for the slow conver- 
gence (resolved radius oc N~^^^) in densi ty profiles of CDM 
haloes simulated at different resolutions iMoore et alJll993 
IChigna et aUEofxl. 

(ix) Most of the affected particles become "relaxed" very 
early and within the first few Gyrs of the evolution of the 
Universe. This is hardest epoch to accurately resolve in a 
cosmological simulation since the relative force errors can 
be large and the densities of forming haloes can be very 
high. 

The high degrees relaxation show that at 2: = many 
particles have completely different energies and orbits com- 
pared to their evolution in the mean field limit (A'^ 00). 
Current cosmological simulations cannot model all the sub- 
tle dynamics like orbital re sonances which can be important 
e.g. during tidal stripping l)Weinberdll99Sf l . But does relax- 
ation also affect the coarse structure of the object, e.g. their 
density profiles? It is unclear as to how one should interpret 
these results for the following reason. The highest rates of 
relaxation are accumulated early during the formation of the 



haloes. Once a subhalo falls into a larger system the particles 
achieve a higher energy and it is not clear that one should 
accumulate the relaxation timescale in the way that we have 
done since the final "hot" system may lose the memory of 
the init ial condition s through violent relaxation processes. 
Indeed. fMoore et alj (.1999') show that the initial conditions 
play little role in determine the final gross structure of dark 
matter haloes. 

The technique of accumulating the relaxation times by 
measuring the local phase density works very well for near 
equilibrium structures. In some instances, the energies of 
particles in a cosmological simulation will increase due to 
the mass increase from merging. From the virial theorem 
the growth in mass must be accompanied by an increase in 
velocity dispersion. Thus it is not clear how relaxation at a 
high redshift propagates to the final time. However we are 
mainly interested in the structure of the central dark mat- 
ter halo and the substructure haloes. These are the regions 
where most effort is going into co mparing theoretical predic- 
tions with observational data ('e.g. lMoorell994l : IStoehr et. all 
l2002h . 

The velocity dispersion in the substructure haloes ac- 
tually decreases slowly with time as mass is stripped. On 
the other hand the velocity dispersion of the central halo 
increases with time as it grows by merging. We can quantify 
the decrease in the mean velocity dispersion in subhaloes or 
the central halo region by extrapolating the particles back- 
wards through time to examine the velocity dispersion in 
the progenitor haloes. For the galaxy cluster, half of the re- 
laxation is accumulated since 2 = 3. We therefore trace all 
the particles in the central few percent of the cluster back 
to 2 = 4 where we find that all the particles lie in the most 
massive collapsed haloes at that time - the most massive 
structures at high redshift form the central cluster region 
through natural biasing. We find that for the central cluster 
particles, a{z — 4)/a{z = 0) = 1.6, i.e. the mean disper- 
sion only increases by 60%. Similarly for cluster substruc- 
ture haloes, the ratio is practically constant through time. 
For this reason we are confident that the accumulated de- 
grees of relaxation are a good indicator of the true relaxation 
of the interesting regions within the simulations. 

We note that the following thought experiment illus- 
trates a possible way in which relaxation in the first haloes 
can affect the central structure of the final halo, even though 
the energies of the particles may have changed over time: 
The first haloes to form are several relaxation times old 
and they achieve this in an near equilibrium state. We have 
demonstrated that our cumulative estimator works very well 
for non-equilibrium systems that have P.E./K.E. ratios as 
large as ten times the equilibrium mean (Figure 5) . We also 
show that these first poorly resolved haloes will develop con- 
stant density cores through relaxation (Figure 8). Now con- 
sider what happens when this halo accretes into a larger 
system. Because of the constant density core the satellite 
will completely disrupt at some distance fr om the centre of 
the final object iMoore. Katz &: Lakell993) . If the accreting 
halo had a steep cusp resolved with more particles then it 
may sink deeper within the po t ential and deposit m ass at 
smaller radn (e.g. 'Barnes"l99gl. ISver fc Whitd ll993). Thus 
the early relaxation could affect the final density profile even 
if the relative energies of the particles were initially quite dif- 
ferent from their final energies. Unfortunately, quantifying 
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this effect and determining if it is at all important is best 
achieved by increasing the resolution by several orders of 
magnitude over the haloes simulated in this paper. 



This paper has been typeset from a TjjX/ I^T^^pC file prepared 
by the author. 



ACKNOWLEDGMENTS 

We thank Adrian Jenkins and Alexander Knebe for helpful 
suggestions and comments. We also would like to thank the 
Swiss supercomputing centre at Manno for computing time 
where many of these numerical simulations were performed. 



REFERENCES 

Barnes J. E., 1999, in Barnes J. E., Sanders D. B., eds, 
Proc. lAU Symp. 186, Galaxy Interactions at Low and 
High Redshift. Kluwer, Dordrecht, p. 137 

Berlin G., 2000, Galactic Dynamics. Cambridge Univ. 
Press, Cambridge 

Bertschinger E., 2001, ApJSS, 137, 1 

Binney J., Knebe A., 2002, MNRAS, 333, 378 

Binney J., Tremaine S., 1987, Galactic Dynamics. Prince- 
ton Univ. Press, Princeton 

Chandrasekhar S., 1942, Principles of Stellar Dynamics. 
Univ. Chicago Press, Chicago 

Fukushige T., Makino J., 2001, ApJ, 557, 533 

Ghigna S., Moore B., Governato F., Lake G, Quinn T., 
Stadel J., 2000, ApJ, 544, 616 

Hayashi E., Navarro J. F., Taylor J. E., Stadel J., Quinn 
T., 2003, ApJ, 584, 541 

Hernquist L., 1990, ApJ, 356, 359 

Huang S., Dubinski J., Carlberg R. G, 1993, ApJ, 404, 73 
Jing Y., Suto Y., 2000, ApJ, 529, L69 

Kazantzidis S., Magorrian J., Moore B., 2004, ApJ, 601, 37 
Klypin A., Kravtsov A. V., Bullock J. S., Primack J. R., 

2001, ApJ, 554, 903 
Knebe A., Green A., Binney J., 2001, MNRAS, 325, 845 
Kravtsov A. V., Klypin A. A., Khokhlov A.M., 1997, ApJS, 

111, 73 

Navarro J. F., Frenk C. S., White S. D. M., 1996, ApJ, 462, 
563 

Moore B., 1994, V.370, Nat, 370, 629 

Moore B., Katz N., Lake G, 1996, ApJ, 457, 455 

Moore B., Governato F., Quinn T., Stadel J., Lake G, 

1998, ApJ, 499, L5 

Moore B., Quinn T., Governato F., Stadel J., Lake G, 

1999, MNRAS, 310, 1147 

Moore B. 2001, In Martel H. & Wheeler J., eds, AIP Conf. 
Proc. 586, The dark matter crisis, p. 73 

Power C, Navarro J. F., Jennkins A., Frenk C. S., White 
S. D. M., 2003, MNRAS, 338, 14 

Quinlan G. D., 1996, New Astronomy, vol. 1, no. 3, 255 

Spitzer L., 1987, Dynamical Evolution of Globular Clus- 
ters. Princeton Series in Astrophysics, Princeton 

Stadel J., 2001, PhD thesis, U. Washington 

Stoehr F., White S. D., Tormen G, Springel V., 2002, MN- 
RAS, 335, L84 

Syer D., White S. D. M., 1998, MNRAS, 293, 337 

Theis C, 1998, A&A, 330, 1180 

Weinberg M. D., 1998, MNRAS, 297, 101 



© 2003 RAS, MNRAS OOO.ITHTTl 



